This is becoming a collection of Python-based example codes focussing on image analysis of biological matter. Apologies for the code often being a little sloppy and "easy-going". Due to a limited time frame, priorities were set to "progress & visualization" instead of "finding the most efficient solution" :)
Whenever a function or method is developed to fulfill the goals, it is outsources into the subfolder "ModulesOwnAndExternalModified", take a look and have fun. Check out my website (www.bio-century.net) for more projects, literature recommendations, links to data bases and further inspirations.
Goal: Count nuclei of a fluorescence micrograph in an automated way
Steps 1: Import micrograph and compute a binary mask.
# sources
# - https://bbbc.broadinstitute.org/BBBC005
# - https://bbbc.broadinstitute.org/search/synthetic%20cells?
# - https://www.geeksforgeeks.org/distance-transformation-in-image-python-opencv/
# - https://www.geeksforgeeks.org/how-to-adjust-the-position-of-a-matplotlib-colorbar/
# - https://opencv24-python-tutorials.readthedocs.io/en/latest/py_tutorials/py_imgproc/py_watershed/py_watershed.html
# - https://scikit-image.org/docs/stable/auto_examples/segmentation/plot_watershed.html
# - https://stackoverflow.com/questions/70361424/extract-white-regions-from-binary-images
# - https://docs.scipy.org/doc/scipy/reference/generated/scipy.ndimage.binary_fill_holes.html
# - https://www.tutorialspoint.com/how-to-compute-the-area-and-perimeter-of-an-image-contour-using-opencv-python
# - https://www.youtube.com/watch?v=WQpXS9gBEu8
# data source
# Synthetic cells, Accession number BBBC005 · Version 1
# A. Lehmussola, P. Ruusuvuori, J. Selinummi, H. Huttunen and O. Yli-Harja,
# "Computational Framework for Simulating Fluorescence Microscope Images With Cell Populations,"
# in IEEE Transactions on Medical Imaging, vol. 26, no. 7, pp. 1010-1016, July 2007, doi: 10.1109/TMI.2007.896925.
#
# A. Lehmussola, P. Ruusuvuori, J. Selinummi, T. Rajala and O. Yli-Harja,
# "Synthetic Images of High-Throughput Microscopy for Validation of Image Analysis Methods,"
# in Proceedings of the IEEE, vol. 96, no. 8, pp. 1348-1360, Aug. 2008, doi: 10.1109/JPROC.2008.925490.
#
# - https://bbbc.broadinstitute.org/BBBC005
# - https://bbbc.broadinstitute.org/search/synthetic%20cells?
# import modules
import cv2
from matplotlib import cm
from matplotlib import pyplot as plt
import matplotlib.pyplot as plt
from mpl_toolkits.axes_grid1 import make_axes_locatable
import numpy as np
from scipy import ndimage
# initialize and define variables and parameter
# fileName = "SIMCEPImages_A23_C96_F1_s19_w1.TIF" # mask
fileName = "SIMCEPImages_A23_C96_F1_s19_1_real.TIF" # "real" image
cellCountTruth = int(fileName[fileName.find('_C')+2:fileName.find('_F')]) # extract real cell number from filename
myThresholdValue = 70 # treshold value for binary mask
myImage = cv2.imread('./Data/' + fileName) # read image
myImageGray = cv2.cvtColor(myImage, cv2.COLOR_BGR2GRAY)
# ret, thresh = cv2.threshold(gray, 0, 1, cv2.THRESH_BINARY+cv2.THRESH_OTSU)
# fill holes
# define binary mask on computed on thresholded image
myImageThresholded = np.zeros(myImageGray.shape)
myImageThresholded[myImageGray > myThresholdValue] = 1
myImageThresholded = ndimage.binary_fill_holes(myImageThresholded).astype(int)
# plot figure
fig, ax = plt.subplots(figsize = (15, 7))
im = ax.imshow(myImageThresholded, cmap = cm.coolwarm)
divider = make_axes_locatable(ax)
cax = divider.new_vertical(size = "5%", pad = 0.3, pack_start = True)
fig.add_axes(cax)
cbar = fig.colorbar(im, cax = cax, ticks = [-1, 0, 1], orientation = 'horizontal')
TitleString = "\"Synthetic cells\": Binary Mask\n"
plt.title(TitleString, fontweight="bold", fontsize = 16, y = 22.3)
plt.suptitle("File: " + fileName, fontsize = 10, y = 0.91)
ax.set_xlabel("binary mask color scheme")
ax.set_xticks([])
ax.set_yticks([])
plt.show()
Step 2: Identify single objects and draw contours
# sources
# - https://www.kaggle.com/code/voglinio/separating-nuclei-masks-using-convexity-defects
# - https://stackoverflow.com/questions/32401806/get-mask-from-contour-with-opencv
# - https://stackoverflow.com/questions/50591442/convert-3-dim-image-array-to-2-dim
# - https://www.tutorialspoint.com/how-to-compute-the-area-and-perimeter-of-an-image-contour-using-opencv-python
# source of main algorithm
# - https://www.kaggle.com/code/voglinio/separating-nuclei-masks-using-convexity-defects
# User: Costas Voglis
# - https://de.mathworks.com/matlabcentral/answers/43435-i-couldn-t-understand-convex-area
# prop.convex_area / prop.filled_area
# convex_area (convex hull): the smallest region that satisfy two conditions: (1) it is convex (2) it contains
# the original region.
# => The ratio increases the less convex the area is (multiple nuclei)
# import modules
import cv2
import numpy as np
import matplotlib.pyplot as plt
from skimage.measure import regionprops
import ModulesOwnAndExternalModified.A_Groundwork as Groundwork
from statistics import mean
# initialize and define variables and parameter
# tresholds for
# 1Cell: single cell
# 2CellCluster: two cells
# gt2CellCluster: greater than two cells
dictThresholdValues = {}
dictThresholdValues["1Cell"] = 0
dictThresholdValues["2CellCluster"] = 1.065
dictThresholdValues["gt2CellCluster"] = 1.195
# compute results
# compute mean area for a single nucleus
meanFilledAreaTmp = Groundwork.meanFilledArea(myImage, myImageThresholded, dictThresholdValues)
# execute main algorith to identify cells and cell cluster
cellCount, contours, myImageContours, myMaskContoursAll = Groundwork.identifyCellContoursAndAreaOverlays(myImage,
myImageThresholded,
dictThresholdValues,
meanFilledAreaTmp)
# print table to compare cellcounts
print("-" * 27, end = '\n')
print(" Number of Cells")
print("-" * 27, end = '\n')
print("Detected: \t\t " + str(len(contours)))
print("True Value: \t\t " + str(cellCountTruth))
print("-" * 27, end = '\n')
print("Merged Cells Delta: \t" + str(len(contours) - cellCountTruth))
# plot figure
fig, ax = plt.subplots(figsize = (15, 7))
im = ax.imshow(myImageContours, cmap = cm.coolwarm)
divider = make_axes_locatable(ax)
cax = divider.new_vertical(size = "5%",
pad = 0.3,
pack_start = True)
TitleString = "\"Synthetic cells\": Determine Cell Contours\n"
plt.title(TitleString, fontweight = "bold", fontsize = 16)
plt.suptitle("File: " + fileName, fontsize = 10, y = 0.91)
ax.set_xticks([])
ax.set_yticks([])
plt.show()
---------------------------
Number of Cells
---------------------------
Detected: 74
True Value: 96
---------------------------
Merged Cells Delta: -22
For illustration purposes: Compute prop.convex_area-to-prop.filled_area distribution on the binary mask, which shall be taken as a criterium for multi-nuclei detection later on and defined in the above section.
# source of main algorithm
# - https://www.kaggle.com/code/voglinio/separating-nuclei-masks-using-convexity-defects
# User: Costas Voglis
# import modules
import cv2
import numpy as np
import matplotlib.pyplot as plt
from skimage.measure import regionprops
from statistics import mean
from matplotlib.gridspec import GridSpec
# contours, hierarchy = cv2.findContours(myImageThresholded.astype(np.uint8), 1, 2)
ratio_conv_filled_All = []
for ii, cnt in enumerate(contours):
M = cv2.moments(cnt)
area = cv2.contourArea(cnt)
perimeter = cv2.arcLength(cnt, True)
perimeter = round(perimeter, 1)
myImageContours = cv2.drawContours(myImage, [cnt], -1, (0, 0, 255), 2)
myMaskZeros = np.zeros(myImage.shape, np.uint8)
myMaskContoursTmp = cv2.drawContours(image = myMaskZeros, contours=[cnt], contourIdx = -1, color = (255, 255, 255), thickness = cv2.FILLED)
myMaskContours = np.zeros(myImage.shape)
myMaskContours[myMaskContoursTmp == 255] = 1
props = regionprops(myMaskContoursTmp, cache = False)
prop = props[0]
ratio_conv_filled_All.append(prop.convex_area / prop.filled_area)
# plot figure
fig = plt.figure(figsize = (16, 5))
gs = GridSpec(nrows = 1, ncols = 1)
ax0 = fig.add_subplot(gs[0, 0])
plt.hist(ratio_conv_filled_All, bins=30)
ax0 = plt.gca()
min_ylim, max_ylim = plt.ylim()
plt.axvline(dictThresholdValues["2CellCluster"], color='k', linestyle='dashed', linewidth = 1)
plt.text(dictThresholdValues["2CellCluster"]*1.005, max_ylim*0.9, 'Borderline 1 or 2 Cells: {:.2f}'.format(dictThresholdValues["2CellCluster"]))
plt.axvline(dictThresholdValues["gt2CellCluster"], color='k', linestyle='dashed', linewidth = 1)
plt.text(dictThresholdValues["gt2CellCluster"]*1.005, max_ylim*0.9, 'Borderline 2 or 3 Cells: {:.2f}'.format(dictThresholdValues["gt2CellCluster"]))
plt.grid(axis = 'y', linestyle = '--', linewidth = 0.3)
ax0.set_title('prop.convex_area / prop.filled_area')
ax0.set_xlabel('prop.convex_area / prop.filled_area')
ax0.set_ylabel('number of nucleic regions')
ax0.set_xticks(np.arange(1, 1.3, 0.01));
ax0.set_yticks(np.arange(0, 30, 1));
Step 3: Identify multi-nuclei regions
# sources
# - https://www.kaggle.com/code/voglinio/separating-nuclei-masks-using-convexity-defects
# - https://stackoverflow.com/questions/32401806/get-mask-from-contour-with-opencv
# - https://stackoverflow.com/questions/50591442/convert-3-dim-image-array-to-2-dim
# - https://www.tutorialspoint.com/how-to-compute-the-area-and-perimeter-of-an-image-contour-using-opencv-python
# compute results
# print table to compare cellcounts
print("-" * 27, end = '\n')
print(" Number of Cells")
print("-" * 27, end = '\n')
print("Detected previously: \t " + str(len(contours)))
print("Detected now: \t\t " + str(cellCount))
print("True Value: \t\t " + str(cellCountTruth))
print("-" * 27, end = '\n')
print("Merged cells delta: \t " + str(cellCount - cellCountTruth))
# plot figure
fig, ax = plt.subplots(figsize = (15, 7))
fig.set_figwidth(12)
fig.set_figheight(12)
ax.imshow(myMaskContoursAll)
TitleString = "\"Synthetic cells\": Detection of Multi-Core Overlays\n"
plt.title(TitleString, fontweight="bold", fontsize = 16)
plt.suptitle("File: " + fileName, fontsize = 10, y = 0.803)
ax.set_xticks([])
ax.set_yticks([])
plt.show()
Clipping input data to the valid range for imshow with RGB data ([0..1] for floats or [0..255] for integers).
---------------------------
Number of Cells
---------------------------
Detected previously: 74
Detected now: 96
True Value: 96
---------------------------
Merged cells delta: 0
Step 4 (add-on): Compute distance transform
# source
# - https://www.geeksforgeeks.org/distance-transformation-in-image-python-opencv/
# import modules
import cv2
import numpy as np
# compute results
distTransform = cv2.distanceTransform(myImageThresholded.astype(np.uint8), cv2.DIST_L2, 3)
fig, ax = plt.subplots(figsize = (15, 7))
im = ax.imshow(distTransform, cmap = cm.coolwarm)
divider = make_axes_locatable(ax)
cax = divider.new_vertical(size = "5%",
pad = 0.3,
pack_start = True)
# plot figure
plt.title("Distance Transform", fontweight="bold", fontsize = 16, y = 1.05)
plt.suptitle("File: " + fileName, fontsize = 10, y = 0.91)
ax.set_xticks([])
ax.set_yticks([])
plt.show()
Step 5 (add-on): Compute watershed results and overlay
# sources
# - https://www.geeksforgeeks.org/distance-transformation-in-image-python-opencv/
# - https://gist.github.com/jgomezdans/402500
# - https://scikit-image.org/docs/stable/auto_examples/segmentation/plot_watershed.html
# - https://stackoverflow.com/questions/31054352/applying-matplotlib-color-map-to-opencv-image
# - https://stackoverflow.com/questions/40895785/using-opencv-to-overlay-transparent-image-onto-another-image
# - https://stackoverflow.com/questions/60765915/python-opencv-floodfill-function-strange-type-errors
# - https://stackoverflow.com/questions/70251361/imextendedmin-and-imimposemin-functions-in-python-for-watershed-seeds-from-dista
# import modules
import codecs, json
import cv2
from matplotlib import cm
import matplotlib.pyplot as plt
import numpy as np
from PIL import Image
from scipy import ndimage
from skimage.segmentation import watershed
# compute results
distTransform_t = ndimage.distance_transform_edt(myImageThresholded.astype(np.uint8))
# plt.imshow(distTransform_t)
waterShedValues = watershed(- distTransform_t)
cv2.imwrite('./Figures/cellscontourized.png', myImageContours)
cv2.imwrite('./Figures/distTransform.png', distTransform * 255)
cv2.imwrite('./Figures/watershedvalues.png', waterShedValues)
image_array = plt.imread("./Figures/watershedvalues.png")
randomValues = np.linspace(0,1,256)
np.random.shuffle(randomValues)
cmap = plt.cm.colors.ListedColormap(plt.cm.jet(randomValues))
dictRandomValues = {}
dictRandomValues["randomValues"] = randomValues.tolist()
with open('./ModulesOwnAndExternalModified/cmap.json', 'w') as myFile:
json.dump(dictRandomValues, myFile)
myFile.close()
imageTmp = Image.fromarray(np.uint8(cmap(image_array)*255))
imageTmp.save('./Figures/watershedvaluesRGB.png')
myBackgroundImage = cv2.imread('./Figures/cellscontourized.png')
myOverlayImage = cv2.imread('./Figures/watershedvaluesRGB.png')
myBackgroundImage = cv2.addWeighted(myBackgroundImage, 0.5, myOverlayImage, 0.5, 0)
# plot figures
fig, axes = plt.subplots(ncols = 2, nrows = 2, figsize = (12, 10), sharex = True, sharey = True)
ax = axes.ravel()
ax[0].imshow(distTransform, cmap = cm.coolwarm)
ax[0].set_title('distance tranform objects')
ax[0].set_xticks([])
ax[0].set_yticks([])
ax[1].imshow(myImageThresholded, cmap = cm.coolwarm)
ax[1].set_title('watershed')
ax[1].set_xticks([])
ax[1].set_yticks([])
ax[2].imshow(myImageContours, cmap = cm.coolwarm)
ax[2].set_title('binary with center of mass')
ax[2].set_xticks([])
ax[2].set_yticks([])
ax[3].imshow(myBackgroundImage, cmap = cm.coolwarm)
ax[3].set_title("Distance Transform of\n" + fileName)
ax[3].set_xticks([])
ax[3].set_yticks([])
plt.show()
Goal: Count nuclei of a set of fluorescence micrographs in an automated way and identify the offset on each image by comparison to the ground truth number
# sources
# - https://www.kaggle.com/code/voglinio/separating-nuclei-masks-using-convexity-defects
# - https://pynative.com/python-list-files-in-a-directory/
# - https://stackoverflow.com/questions/32401806/get-mask-from-contour-with-opencv
# - https://stackoverflow.com/questions/50591442/convert-3-dim-image-array-to-2-dim
# - https://www.tutorialspoint.com/how-to-compute-the-area-and-perimeter-of-an-image-contour-using-opencv-python
# data source
# Synthetic cells, Accession number BBBC005 · Version 1
# A. Lehmussola, P. Ruusuvuori, J. Selinummi, H. Huttunen and O. Yli-Harja,
# "Computational Framework for Simulating Fluorescence Microscope Images With Cell Populations,"
# in IEEE Transactions on Medical Imaging, vol. 26, no. 7, pp. 1010-1016, July 2007, doi: 10.1109/TMI.2007.896925.
#
# A. Lehmussola, P. Ruusuvuori, J. Selinummi, T. Rajala and O. Yli-Harja,
# "Synthetic Images of High-Throughput Microscopy for Validation of Image Analysis Methods,"
# in Proceedings of the IEEE, vol. 96, no. 8, pp. 1348-1360, Aug. 2008, doi: 10.1109/JPROC.2008.925490.
# import modules
import cv2
import numpy as np
import matplotlib.pyplot as plt
import os
from skimage.measure import regionprops
import ModulesOwnAndExternalModified.A_Groundwork as Groundwork
import matplotlib.pyplot as plt
import numpy as np
plotvalue = True
# plotvalue = False
if plotvalue == True:
xpoints = []
ypoints = []
# tresholds for
# 1Cell: single cell
# 2CellCluster: two cells
# gt2CellCluster: greater than two cells
dictThresholdValues = {}
dictThresholdValues["1Cell"] = 0
dictThresholdValues["2CellCluster"] = 1.05
dictThresholdValues["gt2CellCluster"] = 1.2
myThresholdValue = 70
meanFilledAreaAll = []
res = []
dir_path = './DataExternalSyntheticCells'
for path in os.listdir(dir_path):
if os.path.isfile(os.path.join(dir_path, path)):
res.append(path)
# compute mean area for a single nucleus with all file info
for ii, fileNameTmp in enumerate(res):
fileName = fileNameTmp
cellCountTruth = int(fileName[fileName.find('_C') + 2 : fileName.find('_F')])
myImage = cv2.imread('./DataExternalSyntheticCells/' + fileName)
myImageGray = cv2.cvtColor(myImage, cv2.COLOR_BGR2GRAY)
myImageThresholded = np.zeros(myImageGray.shape)
myImageThresholded[myImageGray > myThresholdValue] = 1
myImageThresholded = ndimage.binary_fill_holes(myImageThresholded).astype(int)
meanFilledArea = Groundwork.meanFilledArea(myImage, myImageThresholded, dictThresholdValues)
meanFilledAreaAll.append(meanFilledArea)
meanFilledArea = mean(meanFilledAreaAll)
# execute main algorith on all files to identify cells and cell cluster
print("-" * 138)
print(" " * 20 + "print scheme: " + "<filename> <TrueNumberNuclei> <ComputedNumberNuclei>")
print("-" * 138)
for ii, fileNameTmp in enumerate(res):
fileName = fileNameTmp
cellCountTruth = int(fileName[fileName.find('_C') + 2 : fileName.find('_F')])
myImage = cv2.imread('./DataExternalSyntheticCells/' + fileName)
myImageGray = cv2.cvtColor(myImage, cv2.COLOR_BGR2GRAY)
myImageThresholded = np.zeros(myImageGray.shape)
myImageThresholded[myImageGray > myThresholdValue] = 1
myImageThresholded = ndimage.binary_fill_holes(myImageThresholded).astype(int)
cellCount, contours, myImageContours, myMaskContoursAll = Groundwork.identifyCellContoursAndAreaOverlays(myImage,
myImageThresholded,
dictThresholdValues,
meanFilledArea)
if ii < 6:
print(fileName + ": " + str(cellCountTruth) + " " + str(cellCount) + "\t\t", end = '')
else:
print(fileName + ": " + str(cellCountTruth) + " " + str(cellCount) + "\t", end = '')
if (ii + 1) % 3 == 0:
print("\n", end = '')
xpoints.append(ii + 1)
ypoints.append(cellCountTruth - cellCount)
# plot figure
TitleString = "Quality Check: Automated Nuclei Detection"
x = np.linspace(0, 1, 2)
plt.plot(xpoints, ypoints)
plt.title(TitleString, fontweight="bold", fontsize = 16, y = 1.07)
plt.suptitle("File: " + fileName, fontsize = 10, y = 0.93)
plt.xlabel("File Index")
plt.ylabel("Difference [Truth - Detected_Cellcount]")
plt.xticks(np.arange(0, 72, 5))
plt.yticks(np.arange(0, 3, 1))
plt.show()
------------------------------------------------------------------------------------------------------------------------------------------
print scheme: <filename> <TrueNumberNuclei> <ComputedNumberNuclei>
------------------------------------------------------------------------------------------------------------------------------------------
SIMCEPImages_A01_C1_F1_s01_w1.TIF: 1 1 SIMCEPImages_A01_C1_F1_s02_w1.TIF: 1 1 SIMCEPImages_A01_C1_F1_s03_w1.TIF: 1 1
SIMCEPImages_A02_C5_F1_s01_w1.TIF: 5 5 SIMCEPImages_A02_C5_F1_s02_w1.TIF: 5 5 SIMCEPImages_A02_C5_F1_s03_w1.TIF: 5 5
SIMCEPImages_A03_C10_F1_s01_w1.TIF: 10 10 SIMCEPImages_A03_C10_F1_s02_w1.TIF: 10 10 SIMCEPImages_A03_C10_F1_s03_w1.TIF: 10 9
SIMCEPImages_A04_C14_F1_s01_w1.TIF: 14 14 SIMCEPImages_A04_C14_F1_s02_w1.TIF: 14 14 SIMCEPImages_A04_C14_F1_s03_w1.TIF: 14 14
SIMCEPImages_A05_C18_F1_s01_w1.TIF: 18 18 SIMCEPImages_A05_C18_F1_s02_w1.TIF: 18 18 SIMCEPImages_A05_C18_F1_s03_w1.TIF: 18 18
SIMCEPImages_A06_C23_F1_s01_w1.TIF: 23 23 SIMCEPImages_A06_C23_F1_s02_w1.TIF: 23 23 SIMCEPImages_A06_C23_F1_s03_w1.TIF: 23 23
SIMCEPImages_A07_C27_F1_s01_w1.TIF: 27 27 SIMCEPImages_A07_C27_F1_s02_w1.TIF: 27 27 SIMCEPImages_A07_C27_F1_s03_w1.TIF: 27 27
SIMCEPImages_A08_C31_F1_s01_w1.TIF: 31 31 SIMCEPImages_A08_C31_F1_s02_w1.TIF: 31 31 SIMCEPImages_A08_C31_F1_s03_w1.TIF: 31 31
SIMCEPImages_A09_C35_F1_s02_w1.TIF: 35 35 SIMCEPImages_A09_C35_F1_s03_w1.TIF: 35 35 SIMCEPImages_A10_C40_F1_s02_w1.TIF: 40 40
SIMCEPImages_A10_C40_F1_s03_w1.TIF: 40 40 SIMCEPImages_A11_C44_F1_s01_w1.TIF: 44 44 SIMCEPImages_A11_C44_F1_s02_w1.TIF: 44 44
SIMCEPImages_A11_C44_F1_s03_w1.TIF: 44 44 SIMCEPImages_A12_C48_F1_s01_w1.TIF: 48 48 SIMCEPImages_A12_C48_F1_s02_w1.TIF: 48 48
SIMCEPImages_A12_C48_F1_s03_w1.TIF: 48 48 SIMCEPImages_A13_C53_F1_s01_w1.TIF: 53 53 SIMCEPImages_A13_C53_F1_s03_w1.TIF: 53 53
SIMCEPImages_A14_C57_F1_s01_w1.TIF: 57 57 SIMCEPImages_A14_C57_F1_s02_w1.TIF: 57 57 SIMCEPImages_A14_C57_F1_s03_w1.TIF: 57 57
SIMCEPImages_A15_C61_F1_s01_w1.TIF: 61 61 SIMCEPImages_A15_C61_F1_s02_w1.TIF: 61 61 SIMCEPImages_A15_C61_F1_s03_w1.TIF: 61 61
SIMCEPImages_A16_C66_F1_s01_w1.TIF: 66 65 SIMCEPImages_A16_C66_F1_s02_w1.TIF: 66 66 SIMCEPImages_A16_C66_F1_s03_w1.TIF: 66 66
SIMCEPImages_A17_C70_F1_s01_w1.TIF: 70 70 SIMCEPImages_A17_C70_F1_s02_w1.TIF: 70 69 SIMCEPImages_A17_C70_F1_s03_w1.TIF: 70 70
SIMCEPImages_A18_C74_F1_s02_w1.TIF: 74 74 SIMCEPImages_A18_C74_F1_s03_w1.TIF: 74 73 SIMCEPImages_A19_C78_F1_s01_w1.TIF: 78 78
SIMCEPImages_A19_C78_F1_s02_w1.TIF: 78 78 SIMCEPImages_A19_C78_F1_s03_w1.TIF: 78 78 SIMCEPImages_A20_C83_F1_s01_w1.TIF: 83 83
SIMCEPImages_A20_C83_F1_s03_w1.TIF: 83 82 SIMCEPImages_A21_C87_F1_s01_w1.TIF: 87 85 SIMCEPImages_A21_C87_F1_s02_w1.TIF: 87 87
SIMCEPImages_A21_C87_F1_s03_w1.TIF: 87 87 SIMCEPImages_A22_C91_F1_s02_w1.TIF: 91 91 SIMCEPImages_A22_C91_F1_s03_w1.TIF: 91 91
SIMCEPImages_A23_C96_F1_s01_w1.TIF: 96 95 SIMCEPImages_A23_C96_F1_s02_w1.TIF: 96 95 SIMCEPImages_A23_C96_F1_s03_w1.TIF: 96 95
SIMCEPImages_A24_C100_F1_s01_w1.TIF: 100 100 SIMCEPImages_A24_C100_F1_s02_w1.TIF: 100 100 SIMCEPImages_A24_C100_F1_s03_w1.TIF: 100 99
Goal: Ellipse fit of each nuclei region
# sources
# - https://www.kaggle.com/code/voglinio/separating-nuclei-masks-using-convexity-defects
# - https://de.mathworks.com/matlabcentral/answers/43435-i-couldn-t-understand-convex-area
# - https://stackoverflow.com/questions/32401806/get-mask-from-contour-with-opencv
# - https://stackoverflow.com/questions/50591442/convert-3-dim-image-array-to-2-dim
# - https://stackoverflow.com/questions/62698756/opencv-calculating-orientation-angle-of-major-and-minor-axis-of-ellipse
# - https://www.tutorialspoint.com/how-to-compute-the-area-and-perimeter-of-an-image-contour-using-opencv-python
# - https://www.tutorialspoint.com/how-to-fit-the-ellipse-to-an-object-in-an-image-using-opencv-python
# sources of main algorithms
# https://www.kaggle.com/code/voglinio/separating-nuclei-masks-using-convexity-defects
# User: Costas Voglis
# and
# https://stackoverflow.com/questions/62698756/opencv-calculating-orientation-angle-of-major-and-minor-axis-of-ellipse
# User: fmw42, Fred Weinhaus
# import modules
import cv2
import numpy as np
import matplotlib.pyplot as plt
import ModulesOwnAndExternalModified.A_Groundwork as Groundwork
from skimage.measure import regionprops
import math
fileName = "SIMCEPImages_A23_C96_F1_s02_w1.TIF"
cellCountTruth = int(fileName[fileName.find('_C')+2:fileName.find('_F')])
myThresholdValue = 70
myImage = cv2.imread('./DataExternalSyntheticCells/' + fileName)
myImageGray = cv2.cvtColor(myImage, cv2.COLOR_BGR2GRAY)
myImageThresholded = np.zeros(myImageGray.shape)
myImageThresholded[myImageGray > myThresholdValue] = 1
myImageThresholded = ndimage.binary_fill_holes(myImageThresholded).astype(int)
myImageEllipse = Groundwork.ellipseFit(myImage, myImageThresholded)
# plot figure
fig, ax = plt.subplots(figsize = (20, 10))
im = ax.imshow(myImageEllipse)
divider = make_axes_locatable(ax)
cax = divider.new_vertical(size = "5%",
pad = 0.3,
pack_start = True)
plt.title("Ellipse fit on nuclei", fontweight="bold", fontsize = 16, y = 1.05)
plt.suptitle("File: " + fileName, fontsize = 10, y = 0.91)
ax.set_xticks([])
ax.set_yticks([])
plt.show()
Goal:
Analyse a first real fluorescence dataset and see how well the above algorithms work. This is really to play around with data and get a feeling of the algorithms, while there is no explicit scientific question posed here (due to a lack of data).
Plan:
Import and show data
# sources
# - https://idr.openmicroscopy.org/
# - https://biapol.github.io/blog/robert_haase/browsing_idr/readme.html
# - https://idr.openmicroscopy.org/webclient/api/plates/
# - https://matplotlib.org/stable/gallery/misc/table_demo.html#sphx-glr-gallery-misc-table-demo-py
# - https://matplotlib.org/stable/gallery/subplots_axes_and_figures/subfigures.html
# - https://stackoverflow.com/questions/31726643/how-to-plot-in-multiple-subplots
# - https://towardsdatascience.com/simple-little-tables-with-matplotlib-9780ef5d0bc4
# - https://towardsdatascience.com/plot-organization-in-matplotlib-your-one-stop-guide-if-you-are-reading-this-it-is-probably-f79c2dcbc801
# data source
# https://idr.openmicroscopy.org/webclient/img_detail/3232695/
# 41744_illum_corrected [Well G9, Field 7]
# "https://idr.openmicroscopy.org/webclient/img_detail/3232695/"
# https://idr.openmicroscopy.org/webclient/api/plates/
# "https://idr.openmicroscopy.org/webclient/img_detail/3232695/"
# https://idr.openmicroscopy.org/webclient/api/plates/
# https://idr.openmicroscopy.org/webclient/img_detail/3232695/
# 41744_illum_corrected [Well G9, Field 7]
# "https://idr.openmicroscopy.org/webclient/img_detail/3232695/"
# https://idr.openmicroscopy.org/webclient/api/plates/
# import modules
import json
import matplotlib.pyplot as plt
from mpl_toolkits.axes_grid1 import make_axes_locatable
import numpy as np
import requests
from skimage.io import imread, imshow
import matplotlib.pyplot as plt
from matplotlib.gridspec import GridSpec
import numpy as np
# initialize and define variables and parameter
screen_id = 1751
imageID = "3232695"
myLink = "https://idr.openmicroscopy.org/webclient/render_image/" + imageID + "/"
dictImageMetaData = {}
cell_text = []
# load data
################ Import metadata ################
# configure which dataset to browse
# screen_id = 2302
MAP_URL = "https://idr.openmicroscopy.org/webclient/api/annotations/?type=map&{type}={screen_id}"
# open an interface to the internet
with requests.Session() as session:
# turn the generic MAP_URL into a specific URL for the screen
qs = {'type': 'screen', 'screen_id': screen_id}
qs = {'image_id': imageID}
IMAGE_DETAILS_URL = "https://idr.openmicroscopy.org/webclient/imgData/{image_id}/"
url = IMAGE_DETAILS_URL.format(**qs)
r = session.get(url)
with open("./Data/02_01_ImageMetaData.json", "w") as outfile:
if r.status_code == 200:
json.dump(r.json(), outfile, indent = 4, sort_keys = True)
with open('./Data/02_01_ImageMetaData.json', 'r') as openfile:
dictImageMetaData = json.load(openfile)
# compute results
################ Import image ################
myImage = imread(myLink)
# table 1: file imformation
tableColumnsFileInfo = ["File Info"]
tableRowsFileInfo = ['myLink', 'MAP_URL', 'IMAGE_DETAILS_URL', 'imageId', 'tiles', 'wellId', 'height', 'width']
tableDataFileInfo = [
[str(myLink)],
[str(MAP_URL)],
[str(url)],
[str(dictImageMetaData["meta"]['imageId'])],
[str(dictImageMetaData["tiles"])],
[str(dictImageMetaData["meta"]['wellId'])],
[str(dictImageMetaData["size"]['height'])],
[str(dictImageMetaData["size"]['width'])],
]
cell_text0 = []
for row in range(len(tableDataFileInfo)):
cell_text0.append(tableDataFileInfo[row])
# table 2: channel information
tableColumnsChannelInfo = ['label', 'active status', 'hexacode']
tableRowsChannelInfo = ['Red Channel','Green Channel','Blue Channel','Blue Channel','Blue Channel']
tableDataChannelInfo = [
[str(dictImageMetaData["channels"][0]['label']), str(dictImageMetaData["channels"][0]['active']), str(dictImageMetaData["channels"][0]['color'])],
[str(dictImageMetaData["channels"][1]['label']), str(dictImageMetaData["channels"][1]['active']), str(dictImageMetaData["channels"][1]['color'])],
[str(dictImageMetaData["channels"][2]['label']), str(dictImageMetaData["channels"][2]['active']), str(dictImageMetaData["channels"][2]['color'])],
[str(dictImageMetaData["channels"][3]['label']), str(dictImageMetaData["channels"][3]['active']), str(dictImageMetaData["channels"][3]['color'])],
[str(dictImageMetaData["channels"][4]['label']), str(dictImageMetaData["channels"][4]['active']), str(dictImageMetaData["channels"][4]['color'])],
]
# there is no pixel size available for this micrograph!
cell_text = []
for row in range(len(tableDataChannelInfo)):
cell_text.append(tableDataChannelInfo[row])
rowColors = [
[0.54901961, 0.58559016, 0.77517878, 1. ],
[0.54901961, 0.58559016, 0.77517878, 1. ],
[0.54901961, 0.58559016, 0.77517878, 1. ],
[0.54901961, 0.58559016, 0.77517878, 1. ],
[0.54901961, 0.58559016, 0.77517878, 1. ],
]
colColors = [
[0.54901961, 0.58559016, 0.77517878, 1. ],
[0.54901961, 0.58559016, 0.77517878, 1. ],
[0.54901961, 0.58559016, 0.77517878, 1. ],
]
colColors0 = [
[0.54901961, 0.58559016, 0.77517878, 1. ],
[0.54901961, 0.58559016, 0.77517878, 1. ],
[0.54901961, 0.58559016, 0.77517878, 1. ],
]
rowColors0 = [
[0.54901961, 0.58559016, 0.77517878, 1. ],
[0.54901961, 0.58559016, 0.77517878, 1. ],
[0.54901961, 0.58559016, 0.77517878, 1. ],
[0.54901961, 0.58559016, 0.77517878, 1. ],
[0.54901961, 0.58559016, 0.77517878, 1. ],
[0.54901961, 0.58559016, 0.77517878, 1. ],
[0.54901961, 0.58559016, 0.77517878, 1. ],
[0.54901961, 0.58559016, 0.77517878, 1. ],
]
################ show figure ################
# set up figure grid
fig = plt.figure(figsize = (20, 10))
gs = GridSpec(nrows = 1, ncols = 1)
# plot image
ax0 = fig.add_subplot(gs[0, 0])
ax0.imshow(myImage)
# ax0 = plt.gca()
ax0.get_xaxis().set_visible(False)
ax0.get_yaxis().set_visible(False)# Hide axes border
plt.suptitle("Fluorescence Data", fontweight = "bold", fontsize = 16, y = 0.9)# Add footer
plt.draw()# Create image. plt.savefig ignores figure edge and face colors, so map them.
plt.show()
fig = plt.figure(figsize = (10, 4.5))
gs = GridSpec(nrows = 2, ncols = 1)
ax1 = fig.add_subplot(gs[0, 0])
the_table = plt.table(cellText = cell_text0,
rowLabels = tableRowsFileInfo,
rowColours = rowColors0,
colLabels = tableColumnsFileInfo,
colColours = colColors0,
loc = 'center')
the_table.scale(0.9, 1.4)
ax1 = plt.gca()
ax1.get_xaxis().set_visible(False)
ax1.get_yaxis().set_visible(False)
plt.box(on = None)
ax2 = fig.add_subplot(gs[1, 0])
the_table = plt.table(cellText = cell_text,
rowLabels = tableRowsChannelInfo,
rowColours = rowColors,
colLabels = tableColumnsChannelInfo,
colColours = colColors,
loc = 'center')
the_table.scale(0.9, 1.4)
ax2 = plt.gca()
ax2.get_xaxis().set_visible(False)
ax2.get_yaxis().set_visible(False)
plt.box(on = None)
fig = plt.figure(figsize = (20, 7))
gs = GridSpec(nrows = 1, ncols = 3)
ax3 = fig.add_subplot(gs[0, 0])
plt.imshow(myImage[:,:,0])
ax3 = plt.gca()
ax3.set_title('Hoechst')
ax3.get_xaxis().set_visible(False)
ax3.get_yaxis().set_visible(False)
ax4 = fig.add_subplot(gs[0, 1])
plt.imshow(myImage[:,:,1])
ax4 = plt.gca()
ax4.set_title('ERSyto')
ax4.get_xaxis().set_visible(False)
ax4.get_yaxis().set_visible(False)
ax5 = fig.add_subplot(gs[0, 2])
plt.imshow(myImage[:,:,2])
ax5 = plt.gca()
ax5.set_title('ERSytoBleed')
ax5.get_xaxis().set_visible(False)
ax5.get_yaxis().set_visible(False)
# import modules
import cv2
# save a image using extension
cv2.imwrite("./Data/02_01_myImageRGB.tif", myImage)
cv2.imwrite("./Data/02_01_myImageR.tif", myImage[:,:,0])
cv2.imwrite("./Data/02_01_myImageG.tif", myImage[:,:,1])
cv2.imwrite("./Data/02_01_myImageB.tif", myImage[:,:,2])
True
Define nuclei mask
# import modules
import cv2
from matplotlib import cm
from matplotlib import pyplot as plt
import matplotlib.pyplot as plt
from mpl_toolkits.axes_grid1 import make_axes_locatable
import numpy as np
from scipy import ndimage
# define filenames and parameters
fileName = "02_01_myImageR.tif"
myThresholdValue = 30
myImage = cv2.imread('./Data/' + fileName)
myImageGray = cv2.cvtColor(myImage, cv2.COLOR_BGR2GRAY)
# ret, thresh = cv2.threshold(gray, 0, 1, cv2.THRESH_BINARY + cv2.THRESH_OTSU)
# fill holes
# define binary mask on computed on thresholded image
myImageThresholded = np.zeros(myImageGray.shape)
myImageThresholded[myImageGray > myThresholdValue] = 1
myImageThresholded = ndimage.binary_fill_holes(myImageThresholded).astype(int)
# plot figure
fig = plt.figure(figsize = (20, 10))
gs = GridSpec(nrows = 1, ncols = 3)
ax0 = fig.add_subplot(gs[0, 0])
plt.imshow(myImage[:,:,0])
ax0 = plt.gca()
ax0.set_title('Hoechst')
ax0.get_xaxis().set_visible(False)
ax0.get_yaxis().set_visible(False)
ax1 = fig.add_subplot(gs[0, 1])
plt.imshow(myImageThresholded)
ax1 = plt.gca()
ax1.set_title('Nuclei Binary Mask')
ax1.get_xaxis().set_visible(False)
ax1.get_yaxis().set_visible(False)
Remove small objects from nuclei mask
from skimage import morphology
# sources
# - https://scikit-image.org/docs/dev/api/skimage.morphology.html#skimage.morphology.remove_small_objects
# remove small objects
minimum_pixelarea = 30
myImageThresholdedtmp = []
myImageThresholdedrso = []
myImageThresholdedtmp = myImageThresholded
myImageThresholdedtmp[980+20:1030, 550:585] = 0
myImageThresholdedtmp[980:1025, 550+25:585+20] = 0
myImageThresholdedtmp[990+2:1010, 550:585+20] = 0
myImageThresholdedrso = morphology.remove_small_objects(np.array(myImageThresholdedtmp, dtype=bool), minimum_pixelarea)
myImageRThresholdedRSO = myImageThresholdedrso
# plot figure
fig, ax = plt.subplots(figsize = (15, 7))
im = ax.imshow(myImageThresholdedrso, cmap = cm.coolwarm)
divider = make_axes_locatable(ax)
cax = divider.new_vertical(size = "5%", pad = 0.3, pack_start = True)
fig.add_axes(cax)
cbar = fig.colorbar(im, cax = cax, ticks = [-1, 0, 1], orientation = 'horizontal')
plt.suptitle("Nuclei Binary Mask Cleanup", fontweight = "bold", fontsize = 16, y = 0.92)# Add footer
ax.set_xlabel("binary mask color scheme")
# ax.set_xticks([])
# ax.set_yticks([])
plt.show()
Define mask for cytoplasm and mask out debris
# sources
# - https://scikit-image.org/docs/dev/api/skimage.morphology.html#skimage.morphology.remove_small_objects
# import modules
from skimage import morphology
# define filenames and parameters
fileName = "02_01_myImageB.tif"
myThresholdValue = 10
myImage = cv2.imread('./Data/' + fileName)
myImageGray = cv2.cvtColor(myImage, cv2.COLOR_BGR2GRAY)
# ret, thresh = cv2.threshold(gray, 0, 1, cv2.THRESH_BINARY + cv2.THRESH_OTSU)
# fill holes
# define binary mask on computed on thresholded image
myImageThresholded = np.zeros(myImageGray.shape)
myImageThresholded[myImageGray > myThresholdValue] = 1
myImageThresholded = ndimage.binary_fill_holes(myImageThresholded).astype(int)
# remove small objects
minimum_pixelarea = 30
myImageThresholdedtmp = []
myImageThresholdedrso = []
myImageThresholdedtmp = myImageThresholded
# remove specific debris
myImageThresholdedtmp[980+20:1030, 550:585] = 0
myImageThresholdedtmp[980:1025, 550+25:585 + 20] = 0
myImageThresholdedtmp[990+2:1010, 550:585 + 20] = 0
myImageThresholdedrso = morphology.remove_small_objects(np.array(myImageThresholdedtmp, dtype = bool), minimum_pixelarea)
myImageGThresholdedRSO = myImageThresholdedrso
# plot figure
fig = plt.figure(figsize = (20, 10))
gs = GridSpec(nrows = 1, ncols = 3)
ax0 = fig.add_subplot(gs[0, 0])
plt.imshow(myImage[:,:,0]*2)
ax0 = plt.gca()
ax0.set_title('Cytoplasm')
ax0.get_xaxis().set_visible(False)
ax0.get_yaxis().set_visible(False)
ax1 = fig.add_subplot(gs[0, 1])
plt.imshow(myImageThresholdedrso)
ax1 = plt.gca()
ax1.set_title('Cytoplasm Binary Mask')
ax1.get_xaxis().set_visible(False)
ax1.get_yaxis().set_visible(False)
Review masks
# compute results
contours, myImageContours, myMaskRContoursAll, myMaskRContoursAllLabel = Groundwork.identifyCellContours(myImage, myImageRThresholdedRSO, showCenterOfMass=False, perimeterColor=[20, 240, 20])
contours, myImageContours, myMaskGContoursAll, myMaskGContoursAllLabel = Groundwork.identifyCellContours(myImage, myImageGThresholdedRSO, showCenterOfMass=False, perimeterColor=[20, 240, 20])
contours, myImageContours, myMaskBContoursAll, myMaskBContoursAllLabel = Groundwork.identifyCellContours(myImage, myImageGThresholdedRSO, showCenterOfMass=False, perimeterColor=[20, 20, 240])
myMaskContoursAll = np.zeros(myImage.shape)
myMaskContoursAll[myMaskGContoursAll==1] = 100
myMaskContoursAll[myMaskRContoursAll==1] = 255
print(np.min(myMaskGContoursAllLabel))
print(np.max(myMaskGContoursAllLabel))
# plot figure
fig = plt.figure(figsize = (20, 19))
gs = GridSpec(nrows = 1, ncols = 3)
ax0 = fig.add_subplot(gs[0, 0])
plt.imshow(myImageContours*2)
ax0 = plt.gca()
ax0.set_title('Cells With Boundaries')
ax0.get_xaxis().set_visible(False)
ax0.get_yaxis().set_visible(False)
ax1 = fig.add_subplot(gs[0, 1])
plt.imshow(myMaskContoursAll.astype(np.uint8))
ax1 = plt.gca()
ax1.set_title('Cell Binary Masks')
ax1.get_xaxis().set_visible(False)
ax1.get_yaxis().set_visible(False)
ax2 = fig.add_subplot(gs[0, 2])
plt.imshow(myMaskGContoursAllLabel.astype(np.uint8))
ax2 = plt.gca()
ax2.set_title('Cell Binary Masks with Labels')
ax2.get_xaxis().set_visible(False)
ax2.get_yaxis().set_visible(False)
0.0 78.0
Eliminate everything without a nucleus from cytoplasmic mask
# import modules
import math
# define filenames and parameters
plotstatus = False
jj=0
myMaskContoursAllCorrected = np.zeros((myImage.shape[0], myImage.shape[1],int(np.max(myMaskBContoursAllLabel))))
print(myMaskContoursAllCorrected.shape)
# compute results
for ii in range(1,int(np.max(myMaskBContoursAllLabel))+1):
myMaskTmp = np.zeros((myImage.shape[0], myImage.shape[1]))
myMaskTmp[myMaskBContoursAllLabel==ii]=1
if np.max(myImageRThresholdedRSO + myMaskTmp) > 1:
print(str(ii) + ", ", end='')
myMaskContoursAllCorrected[:,:,ii-1] = myMaskTmp
# plot figure
if plotstatus==True:
f, axarr = plt.subplots(20,5,figsize = (20, 100))
for ii in range(int(np.max(myMaskBContoursAllLabel))):
xindex = ii % 5
yindex = math.floor(ii / 5)
axarr[yindex,xindex].imshow(myMaskContoursAllCorrected[:,:,ii])
axarr[yindex,xindex].set_title(ii)
(1080, 1080, 78) 3, 4, 6, 7, 8, 9, 10, 11, 12, 14, 15, 17, 18, 20, 22, 23, 24, 25, 26, 27, 28, 29, 30, 31, 32, 34, 35, 36, 37, 39, 41, 42, 43, 44, 45, 46, 47, 49, 50, 52, 53, 54, 55, 56, 57, 58, 59, 61, 62, 63, 64, 65, 66, 67, 68, 69, 70, 71, 76, 77, 78,
Review masks
# compute results
result = np.sum(myMaskContoursAllCorrected[:,:,1:], axis = 2)
myMaskContoursAll_Marks = cv2.imread('./Figures/myMaskContoursAll_Marks.png')
# plot figure
fig = plt.figure(figsize = (16, 15))
gs = GridSpec(nrows = 2, ncols = 2)
ax0 = fig.add_subplot(gs[0, 0])
plt.imshow(myMaskContoursAll_Marks.astype(np.uint8))
ax0 = plt.gca()
ax0.set_title('Cell Binary Masks')
ax0.get_xaxis().set_visible(False)
ax0.get_yaxis().set_visible(False)
cv2.imwrite('./Figures/myMaskContoursAll.png', myMaskContoursAll)
ax1 = fig.add_subplot(gs[0, 1])
plt.imshow(result)
ax1 = plt.gca()
ax1.set_title('Cell Binary Masks after Cleanup')
ax1.get_xaxis().set_visible(False)
ax1.get_yaxis().set_visible(False)
Ellipse fit on nuclei and whole cell
# sources
# - https://www.kaggle.com/code/voglinio/separating-nuclei-masks-using-convexity-defects
# - https://de.mathworks.com/matlabcentral/answers/43435-i-couldn-t-understand-convex-area
# - https://stackoverflow.com/questions/32401806/get-mask-from-contour-with-opencv
# - https://stackoverflow.com/questions/50591442/convert-3-dim-image-array-to-2-dim
# - https://stackoverflow.com/questions/62698756/opencv-calculating-orientation-angle-of-major-and-minor-axis-of-ellipse
# - https://www.tutorialspoint.com/how-to-compute-the-area-and-perimeter-of-an-image-contour-using-opencv-python
# - https://www.tutorialspoint.com/how-to-fit-the-ellipse-to-an-object-in-an-image-using-opencv-python
# sources of main algorithms
# https://www.kaggle.com/code/voglinio/separating-nuclei-masks-using-convexity-defects
# User: Costas Voglis
# and
# https://stackoverflow.com/questions/62698756/opencv-calculating-orientation-angle-of-major-and-minor-axis-of-ellipse
# User: fmw42, Fred Weinhaus
# import modules
import cv2
import numpy as np
import matplotlib.pyplot as plt
import ModulesOwnAndExternalModified.A_Groundwork as Groundwork
from skimage.measure import regionprops
import math
fileName = "02_01_myImageRGB.tif"
myImage = cv2.imread('./Data/' + fileName)
# compute results
myImageEllipseR = Groundwork.ellipseFit(myMaskContoursAll.astype(np.uint8), myImageRThresholdedRSO)
myImageEllipse = Groundwork.ellipseFit(myMaskContoursAll.astype(np.uint8), result)
# plot figure
fig = plt.figure(figsize = (16, 15))
gs = GridSpec(nrows = 1, ncols = 2)
ax0 = fig.add_subplot(gs[0, 0])
plt.imshow(myImageEllipseR)
ax0 = plt.gca()
ax0.set_title('Ellipse Fit on Nucei')
ax0.get_xaxis().set_visible(False)
ax0.get_yaxis().set_visible(False)
ax1 = fig.add_subplot(gs[0, 1])
plt.imshow(myImageEllipse)
ax1 = plt.gca()
ax1.set_title('Ellipse Fit on Whole Cells')
ax1.get_xaxis().set_visible(False)
ax1.get_yaxis().set_visible(False)
Compute Mean Area of Nuclei and Mean Intensity
# import modules
import math
# define filenames and parameters
plotstatus = False
jj = 0
myMaskRContoursAllCorrected = np.zeros((myImage.shape[0], myImage.shape[1],int(np.max(myMaskRContoursAllLabel))))
myMaskTmpIntensity = np.zeros((myImage.shape[0], myImage.shape[1],int(np.max(myMaskRContoursAllLabel))))
filledAreaNuclei = []
meanIntensityNuclei = []
print(myMaskContoursAllCorrected.shape)
# compute results
for ii in range(1,int(np.max(myMaskRContoursAllLabel))+1):
myMaskTmp = np.zeros((myImage.shape[0], myImage.shape[1]))
myMaskTmp[myMaskRContoursAllLabel==ii]=1
if np.max(myImageRThresholdedRSO + myMaskTmp) > 1:
print(str(ii) + ", ", end = '')
myMaskRContoursAllCorrected[:,:,ii-1] = myMaskTmp
numberOfPixel = np.sum(np.sum(myMaskTmp, axis = 0), axis = 0)
filledAreaNuclei.append(numberOfPixel)
myMaskTmpIntensity[:,:,ii-1]=myMaskTmp*myImage[:,:,2]
integratedIntensity = np.sum(np.sum(myMaskTmp*myImage[:,:,2], axis = 0), axis = 0)
meanIntensityNuclei.append(integratedIntensity / numberOfPixel)
# plot figure
if plotstatus==True:
f, axarr = plt.subplots(17,5,figsize = (20, 100))
for ii in range(int(np.max(myMaskRContoursAllLabel))):
xindex = ii % 5
yindex = math.floor(ii / 5)
axarr[yindex,xindex].imshow(myMaskTmpIntensity[:,:,ii])
axarr[yindex,xindex].set_title(ii)
print("")
print("filledAreaNuclei", filledAreaNuclei)
print("meanIntensityNuclei", meanIntensityNuclei)
(1080, 1080, 78) 1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11, 12, 13, 14, 15, 16, 17, 18, 19, 20, 21, 22, 23, 24, 25, 26, 27, 28, 29, 30, 31, 32, 33, 34, 35, 36, 37, 38, 39, 40, 41, 42, 43, 44, 45, 46, 47, 48, 49, 50, 51, 52, 53, 54, 55, 56, 57, 58, 59, 60, 61, 62, 63, 64, 65, 66, 67, 68, 69, 70, 71, 72, 73, 74, 75, 76, 77, 78, 79, 80, 81, 82, 83, 84, filledAreaNuclei [501.0, 683.0, 1072.0, 936.0, 988.0, 913.0, 739.0, 778.0, 1035.0, 1056.0, 666.0, 705.0, 713.0, 838.0, 1016.0, 927.0, 1020.0, 886.0, 747.0, 862.0, 801.0, 1012.0, 770.0, 1169.0, 775.0, 775.0, 1125.0, 714.0, 786.0, 247.0, 873.0, 1059.0, 1230.0, 1076.0, 967.0, 749.0, 669.0, 954.0, 1002.0, 883.0, 664.0, 30.0, 769.0, 835.0, 1134.0, 914.0, 907.0, 671.0, 757.0, 592.0, 678.0, 868.0, 780.0, 777.0, 636.0, 641.0, 873.0, 341.0, 767.0, 648.0, 706.0, 629.0, 470.0, 843.0, 883.0, 1062.0, 608.0, 733.0, 1447.0, 386.0, 243.0, 766.0, 1103.0, 914.0, 967.0, 110.0, 910.0, 1105.0, 638.0, 812.0, 1076.0, 614.0, 173.0, 390.0] meanIntensityNuclei [27.744510978043913, 31.275256222547583, 30.090485074626866, 35.36645299145299, 31.772267206477732, 30.428258488499452, 25.882273342354534, 29.443444730077122, 33.19420289855073, 30.334280303030305, 31.06756756756757, 31.78581560283688, 33.92145862552594, 33.582338902147974, 36.56496062992126, 47.49514563106796, 35.19803921568627, 32.20428893905192, 30.97590361445783, 29.858468677494198, 35.90761548064919, 41.30335968379447, 27.41948051948052, 35.221556886227546, 28.032258064516128, 29.52516129032258, 31.232, 28.84173669467787, 29.0381679389313, 25.453441295546558, 31.079037800687285, 32.271010387157695, 29.70081300813008, 37.69423791821561, 37.26059979317477, 35.30574098798398, 32.18236173393124, 30.947589098532493, 36.46906187624751, 36.20271800679502, 35.42168674698795, 24.9, 27.20416124837451, 33.41676646706587, 35.62522045855379, 33.26258205689278, 42.55678059536935, 38.80923994038748, 27.21664464993395, 31.33783783783784, 30.781710914454276, 33.02304147465438, 35.87307692307692, 24.82110682110682, 29.663522012578618, 33.037441497659906, 27.54982817869416, 24.35777126099707, 34.333767926988266, 29.391975308641975, 33.446175637393765, 32.57551669316375, 30.35531914893617, 26.018979833926455, 26.13929784824462, 31.45762711864407, 35.73848684210526, 33.83765347885402, 35.57567380787837, 33.88341968911917, 28.02880658436214, 30.97911227154047, 40.07978241160471, 44.78774617067834, 36.847983453981385, 33.518181818181816, 35.46263736263736, 31.249773755656108, 33.06426332288401, 55.491379310344826, 38.187732342007436, 27.64332247557003, 117.09826589595376, 86.1]
Plot Mean Area of Nuclei and Mean Intensity
# plot figure
fig = plt.figure(figsize = (16, 5))
gs = GridSpec(nrows = 1, ncols = 2)
ax0 = fig.add_subplot(gs[0, 0])
plt.hist(filledAreaNuclei, bins=30)
ax0 = plt.gca()
ax0.set_title('Mean Area Distribution of Nuclei')
ax0.set_xlabel('mean area / px')
ax0.set_ylabel('number of nuclei')
ax1 = fig.add_subplot(gs[0, 1])
plt.hist(meanIntensityNuclei, bins=30)
ax1 = plt.gca()
ax1.set_title('Mean Intensity Distribution of Nuclei')
ax1.set_xlabel('mean intensity / au')
ax1.set_ylabel('number of nuclei')
Text(0, 0.5, 'number of nuclei')
Goal: Segmentation (first approach) using the Voronoi tesselation.
Steps:
masterarray, myImageInclCenters = Groundwork.ellipseFitCenterPos(myMaskContoursAll.astype(np.uint8), myImageRThresholdedRSO)
# print(masterarray)
import sys
import numpy as np
from scipy.spatial import Voronoi, voronoi_plot_2d
from matplotlib import pyplot as plt
from skimage import io, filters
from skimage.feature import peak_local_max
fileName = "02_01_myImageRGB.tif"
myImage = cv2.imread('./Data/' + fileName)
fig = plt.figure(figsize = (15, 10))
ax = fig.add_subplot(111)
ax.imshow(ndimage.rotate(myImageInclCenters, 0))
cv2.circle(myImage, (int(masterarray[3, 0]), int(masterarray[3, 1])), 5, (255, 0, 0), -1)
vor = Voronoi(masterarray)
ax.set_xticks([])
ax.set_yticks([])
plt.show()
# sources
# - https://stackoverflow.com/questions/20515554/colorize-voronoi-diagram
# - https://stackoverflow.com/questions/52184271/image-and-voronoi-diagram-on-the-same-figure
# import modules
import matplotlib.pyplot as plt
import ModulesOwnAndExternalModified.A_Groundwork as Groundwork
import numpy as np
from scipy.spatial import Voronoi
# compute Voronoi tesselation
vor = Voronoi(masterarray)
# plot
regions, vertices = Groundwork.voronoi_finite_polygons_2d(vor)
fig = plt.figure(figsize=(15, 10))
ax = fig.add_subplot(111)
for region in regions:
polygon = vertices[region]
plt.fill(*zip(*polygon), alpha=0.4)
plt.plot(masterarray[:,0], masterarray[:,1], 'ko')
ax.set_xticks([])
ax.set_yticks([])
plt.imshow(ndimage.rotate(myImageInclCenters, 0))
plt.show()
The results on the given dataset are ok but not remarkable. Still, this algorithm is very robust and one can work with this segmentation: Each of the ellipse centers serves as a seed and the entire set is radially expanded. The points of contact with the neighboring areas limit the evaluation area for each cell. The real strength of the approach is when it comes to dense (i.e. confluent) cell or tissue recordings, where there is hardly any background signal.
# export code to html
import os
os.system("jupyter nbconvert AnalysisOfBioImages.ipynb --to html")
0